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ABSTRACT 


This report treats the problem of mathematically defining a 
smooth surface, z = f(x,y), passing through a finite set of given points 
(xi,yi,Zi, i = In particular, it is not assumed that the 

given (x^,yj^) values lie in any special pattern such as at the nodes 
of a rectangular grid. 

The literature relating to this problem is briefly reviewed. 
An algorithm is described that first constructs a triangular grid in 
the (x,y) domain, next estimates first partial derivatives at the nodal 
points, and finally does interpolation in the triangular cells using 
a method that gives continuity overall. Performance of software 
implementing this algorithm is discussed. 

New theoretical results are presented that provide valuable 
guidance in the development of algorithms for constructing triangular 
grids. 


IV 
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SECTION 1 

PHOPLEM STATEMENT AND GUIDE TO THE LITER ATUhE 


1 . 1 INTRODUCTION 

This paper is a r-esult of our fourth effort in software 
for surface representation. We developed subroutines for rectangular 
grid contour plotting in 1965 with N. Block and R. Garrett, least squares 
bicubic spline surface fitting in 1970 with R. hanson, and contour 
plotting via triangular grid construction and linear interpolation 
in 1972. 


The letter two subroutines deal with irregularly located 
data, however, applications continue to arise in which one would like 
the interpolatory capability of the triangular grid program but with 
at least C^ continuity. Such an algorithm with underlying theory and 
implementing software are the topics of this paper. 

In Section I we introduce the problem and give a brief 
survey of the pertinent literature. Section II describes our algorithm 
and concludes with examples of surfaces produced by our new subroutines. 
We express appreciation to Bob Barnhill and Frank Little for valuable 
discussions that particularly influenced our triangulation algorithm 
of Section 2.2. 

There has been practically no theory to guide the development 
of algorithms for triangulation and no practical static global criterion 
to characterize a preferred triangulation. We are indebted to Michael 
Powell and Robin Sibson for conversations and correspondence in 1976 
that introduced us to Thiessen proximity regions and the fact that 
this concept can be used to define a triangulation as is related in 
Section 3.1.2. 

In our initial effort to determine the relationship of 
the Thiessen criterion to the max-min angle criterion we had used in 
1972, we discovered the circle criterion, which served as a convenient 
mathematical link oetween the other two. The outcome is the material 
of Section 3.1, showing the equivalence of these three criteria when 
used for local optimization of a triangular grid. 

The local equivalence results opened the way to certain 
global equivalences reported in Sec* ion 3.2 and new algorithmic insights 
reported in Sections 3.3 and 3.^*. 

Our conclusions regarding the state of the art for this 
problem appear in Section 3.5. 


1.2 PROBLEM STATEMENT 

The following surface interpolation problem will be treated: 
Given a set of triples of data (xj, y^ , z^), i s 1, .... n, construct 
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a conveniently computable Cl function f(x,y) satisfying the Interpolation 
conditions 


Zi s f(Xi, y^), i s 1,...,n 


The data (xj^, y^) are not assumed to lie in any special pattern 
as at the nodes of a rectangular grid. It is assumed* however, 
all (xi, yi) pairs are distinct; i.e., (xj., y^) s (xj, yj) only 
i = j • 


such 

that 

if 


1.3 EXPECTED APPLICATIONS 

The usual situation in which the author has seen a need 
for this type of computation is that in which a scientist or engineer 
has in hand a set of (xj, y^, data representing measured or computed 
values of some phenomenon and desires to obtain a visual impression 
of a smooth surface of the form z = f(x, y) interpolating the data. 

In such a case, an Interpolation algorithm, such as is treated in this 
paper, must be interfaced with algorithms for contour plotting or surface 
perspective plotting. If, as is the case at JPL, subroutines are available 
for doing contour or surface perspective plotting for data given on a 
rectangular grid, then the surface interpolation algorithm can be used 
to produce the values needed at the lattice points of a rectangular grid. 

Other applications have arisen which can be regarded as 
the inverse of contour plotting. Certain handbook data is available 
in the form of contour plots. To use the data in a computer program 
it is necessary to produce a computable representation of the function 
depicted by the contour plots. A convenient way to do this is to develop 
a list of (xj^, yj^, values from appropriately spaced points along 
the contour lines and then use a surface Interpolation algorithm such 
as is discussed in this paper. 

We have also seen applications which can be regarded as 
implicit function problems. One may have a rectangular table or a 
contour plot giving z as a function of x and y, but then need to be 
able to determine x as a function of y and z in some computational 
procedure. If the data has appropriate monotonicity for this to make 
sense, then the interpolation algorithm of this paper can be used foi 
such problems. 


1.4 PUBLISHED WORK ON SURFACE INTERPOLATION TO IRREGULARLY 

LOCATED DATA 

A variety of algorithmic ideas have been developed for 
this problem or closely related problems. 

Two of the most recent papers giving methods for C1 surface 
interpolation to irregularly located data are Akima (1975) and McLain 
(1976). Akima ’s report contains listings of a set of Fortran subroutines 


1-2 



77-30 


to handle this problem. Ihis cods and a second version of it using 
a more economical triangulation subroutine due to Lawson (1972) have 
been made available to requestors by Akima. 

both Akima (1975) and McLain (1976) contain introductory 
sections giving useful brief characterizations of other approaches, 
particularly those of Bengtsson and Nordbeck (1964), Shepard (1968), 
Maude (1973), and McLain (1974). Franke (1975) reports on computer 
tests of eleven Jiethods constructed using a combination of ideas from 
Sard (1963), Mansfielc. (1972), Maude (1973), McLain (1974), Nielson 
(1974), and Barnhill and Nielson (1974). 

Powell (1976) and ochuiraker (1976) give surveys of methods 
for surface fitting and related pn-blems of bivariate function repre- 
sentation . 


The computerized representation of surfaces is a central 
issue in the field of computer-aided geometric design (CAGD) and plays 
an important role in the field ot finite element methods (FEM). For 
discussions of surface repreoentation from the point of view of CAGD 
see Forrest (1972) and Barnhill (1977). For descriptions of surface 
elements used in FEM see Birkhoff and Mansfield (1974) as well as any 
of the current books on FEM. 

Some methods of building a triangular grid with a given 
set of nodes start by locating the boundary points of the convex :<ull 
of the point set. Algorithms for locating these boundary points are 
treated in Graham (1972), Jarvis (1973), Preparata and Hong (1977), 
and Eddy (1976). 

Some interesting new triangular grid elements providing 
continuity through the use of piecewise quadratics are described 
in Powell and Sabin (1976). 
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SECTION II 

AN ALORITHM FOR TRIANGULATION AND SURFACE INTERPOLATION 


2.1 OUTLINE OF IhE ALGORITHMIC APPROACH SELECTED 

Our approach to the C^ surface interpolation problem consists 
of the following four steps. 

1 ) Construct a triangular gria covering the convex hull 
of the given set of (xj, yj) data using the (xj^, y^) 
data points as vertices of the triangular cells. 

(2) Estimate first partial derivatives of z with respect 
to X and y at each of the (xj, y^) data points. 

(3) For an arbitrary (x, y) point, perform a lookup in 
the triangular grid to identify which triangle, if 
any, contains the point. 

(4) For an arbitrary (x, y) point in a triangle, compute 
an interpolated value of z and optionally of dz/dx 
and 3z/dy also. Make use of the nine items of data 
associated with the triangle, i.e., the values of 

Zi and its two first partial derivatives at each 
of the three vertices. 

This top level description of the approach also characterizes 
the methods of Akima (1975) and McLain (1976), with the exception that 
their methods estimate different quantities at Step 2 for use in the 
interpolation at Step 4. At the more detailed level of devising algorithms 
for each of the four steps, there are substantial differences between 
our approach and that of Akima and that of McLain. 

The four steps will be discussed in the following four 

sections. 


2.2 CONSTRUCTING A TRIANGULAR GRID 

Given the set S of distinct points, (xj^, y^), i : 1,..., n, 
the triangular grid T to be constructed is to cover the convex hull 
of this set of points. Each triangular cell In the grid Is to have 
three of the given points as Its vertices and Is to contain no other 
points of S as Interior or boundary points. 

Conceptually there Is no difficulty In manually constructing 
such a triangular grid. For example, one can Just start, drawing edges 
connecting pairs of points and continue as long as edges can be drawn 
that do not intersect any previously drawn edges. An edge must not 
contain any points of S other than its own two endpoints. 
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In general, there exist many different tri angulations r.i 
a set S. It is noteworthy, however, that all possible triangulatior^ 
of S have the same number of triangles and the same number of edges 
Let nt) denote the number of points of S on the boundary of the convex 
hull of S and let n^ denote the number in the interior so that 
n s n(} ♦ n^> Then the number of triangles is 


n^ s nfc ♦ 2 ( ni - 1 ) 1 2n 


and the number of edges is 


n^ r 2nb ♦ 3 (n^ - 1) 3n 


Taking these relationships into account we selected the 
data structure illustrated by Figures 2-1 and 2-2 to represent a tri- 
angular grid. Column 1 of Figure 2-1 is not stored so each triangle 
is represented in storage by six integers. Since n^ 1 2n, a total 
of 12n integer storage locations suffices to represent the triangular 
grid for n points. This of course is in addition to storage for the 
(xi, yi) data. 

Three sjbroutlnes for storing and fetching in this structure 
are used throughout the package so that the actual mode of storing 
these pointers can be 'hidden'' from the main subroutines. On many 
computers one could easily pack two or three of these pointers per 
computer word Our triangulation algorithm has the property that it 
makes additions and changes to the list of triangles but no deletions, 
'ihus there is no garbage collection problem. 


TRIANGLE 

INDEX 

INDICES OF ADJACENT TRIANGLES 
IN COUNTERCLOCKWISE ORDER. 

A ZERO INDICATES THE REGION 
EXTERIOR TO THE TRIANGULAR 
GRID 

INDICES OF VERTEX POINTS 
IN COUNTERCLOCKWISE 
ORDER. THE F:'^ '>T VERTEX IS 
AT THE POINT OF CONTACT 
OF THE THIRD AND FIRST 
NEIGHBORING TRIANGLE 

1 

2 

0 

4 

5 

8 

7 

2 

5 

3 

1 

5 

3 

8 

3 

• 
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• 

0 

« 
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3 

• 
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• 
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• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

* 

• 
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Figure 2-1. Data Structure hepresenting a Triangular Grid 
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Figure 2-2. The Portion of a Triangular 
Grid Described by the Data 
Structure of Fig. 2-1 


In some reasonable triangulation algorithms the number 
of boundary points of the convex hulls of a sequence of subsets of 
S has an effect on the operation count. Thus it is desirable to have 
some notion of the expected value of n^ as a function of n. Clearly, 
any value of n^ from 3 to n is possible. 

Consider the case in which S is a random sample of n points 
from a uniform distribution on a disc. Let v^j denote the expected 
value of nb in this case. Renyi and Sulanki (1963, 1964) and also 
Rayna>’'1 (197C; show that Vfj is asymtotically proportional to n^^S as 
r — OD. Efron (1965) derives the following formula for v^: 
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He have evaluated this formula by numerical Integration 
and corroborated the results by a cranputer program that counted the 
number of boundary points of the convex hulls of a large number of 
point sets generated as pseudorandom samples from a uniform distribution 
on the unit disc. Selected values are shown in Table 2>1. 

Since a given point set S generally admits many different 
triangulations, how is one to characterize and produce a preferred 
triangulation? So far there does not appear to be any fundamentally 
compelling best answer to this question. 

A very satisfactory candidate, however, emerges from the 
theoretical results presented in Section III. There it is shown that 
three differently stated criteria i ^r choosing the preferred triangulation 
cf a quadrilateral are in fact equivalent in that they produce the 
same decisioas in all cases. It is also shown that if all quadrilaterals 
r«nslsting of pairs of adjacent triangles in a triangulation satisfy 
one of these optimality criteria, then the triangulation as a »*hole 
has some pleasant properties, and in fact is unique to within some 
arbitrariness associated with subsets of four or more neighboring points 
lying on a common circle. 

It is further shown that these local criteria have favorable 
properties for use in triangulation algorithms. In particular it is 
shown that the final trlangulatlon is reached in a finite number of 
steps and that the operation of changing an optimized triangulaticn 
to include a new data point has some properties that can be exploited 
to simplify an algorithm. 

Our triangulation subroutine TPIGRD presently uses the 
max-min angle criterion as its local optimization procedure. This 
is one of the three equivalent local criteria defined in Section 3.1. 


Table 2-1. Vjj Is the Expected Number of Boundary Points in the 
Convex Hull of an n-Point Sample from the Uniform 
Distribution on a Disc 


n 

Vn 

cn = Vn/n1/^ 

4 

3.7 

2.3 

10 

6.1 

2.8 

100 

15.2 

3.3 

1000 

33.6 

3.4 

10000 

72.6 

3.4 
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The TRIGRD algorithm starts by finding the point of S having 
the smallest x coordinate. Any ties are broken by minimizing on y. 

The point p* found in this way is an extreme point of the convex hull 
of S. Finding p* requires 0 (n) operations. 

The points of S are next sorted on increasing (squared) 
Euclidean distance from p*. Denote the points with this ordering by 
q^, q 2 » .... q^ with q^ = p*. We estimate the operation count for 
the sort to be 0 (n log n). 

The first edge is drawn connecting q^ and q 2 . The next 
point in the q^ sequence not colinear with qi and q 2 is connected to 
qi and q 2 to form the first triangle. If this third vertex is not 
q 3 but rather qj^ with k > 3, relabel the points q 3 through qi{ so that 
q)( is called q 3 and the indices of the intervening points are each 
increased by one. These steps assure that qj is strictly outside the 
convex hull of ^q^. .... for all j * H. 3 . .... 

Let c denote the centroid of Aqiq 2 q 3 . Let r be the half 
ray from c through q^. When an angular coordinate is needed for a 
point, the angle will be measured counterclockwise around c from r. 

Note that the angular coordinate of q^ is zero, and all other points 
q^ for i > 1 have angular coordinates strictly between 0 and 2ir. The 
program does not actually compute angles but rather computes a less 
expensive function monotonlcally related to the angle. 


Build an initial boundary list consisting of q^. q 2 . Q 3 . 
and qi (again) along with their angles, assigning the angle 2 ir to the 
second occurrence of qi . 


This finishes the preliminaries. The algorithm proceeds to loop 
through the points qjj, k = 4, ..., n, doing the following for each one: 


Determine the angular coordinate of qi^ and use that coordinate 
as a key to search for a pair of boundary points whose angles bracket 
that angle. This is a linear search requiring an average of nj,^*^V2 
scalar comparisons, where n^^^) is the number of points on the boundary 
of the convex hull of |qi, .... qic_if. If we estimate to be 

about 3 k ^^3 (recall Table 2 - 1 ). then the total cost of this lookup 
as k runs from 4 to n is 0 (n**^ 3 ). This appears to be the highest-order 
operation count in the triangulation algorithm. 


Having found two boundary points to which q|( can be connected, 
attach q)( to these points and record the edge opposite q)( in the new 
triangle in a stack, identifying edges to be tested for possible swapping. 

If the stack is nonempty, unstack one edge and apply the 
local optimization procedure to it. If the decision is to swap the 
edge, do so, and stack the two edges that are opposite q)( in the two 
new triangles. Continue processing the stack until it is empty. 


When the stack is empty try to connect q|( to another neighboring 
boundary point. If this is possible, then run through the stacking 
and testing again, starting with the edge opposite in the new triangle. 
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When cannot be connected to any more boundary points, the processing 
of q|( Is completed. 

We estimate the average operation count to process q^^ is 
a constant. Independent of k. Thus the total cost to process all the 
points q|j, k = 4, . . . , n Is 0(n). 

The total operation count for TRIGRD Is thus estimated 
to be 0(n^/3) + o(n log u) ♦ 0(n). Actual timing on the Uni vac 1108 
for cases having four different values of n are shown in Table 2-2. 

The data of Table 2-2 suggests that in the range 25 < n < 500, 
either OCn^'^^) or 0(n log n) may be used as a model of the execution 
time of this trlangulatlon algorithm. 


2.3 ESTIMATING PARTIAL DERIVATIVES AT THE GRID NODES 

To estimate 9z/9x and 9z/9y at a nodal point p = (xj^, yj^) 
of the triangular grid, the subroutine ESTPD sets up and solves a local 
least squares problem. All of the immediate grid neighbors of point 
p are used up to a maximum of 16 immediate neighbors. If the number 
of immediate neighbors is less than six, then additional nodes beyond 
the immediate neighbors are used to bring the total number of points 
in addition to p up to six. 

A six-parameter quadratic polynominal in x and y is fit 
to the z data values at this set of points. The quadratic is forced 
to interpolate the z value at p, _nd it fits the remaining points in 
a weight^-J least squares sense. The weighting is used to diminish 
the effect f the more distant points. 


Table 2-2. t Denotes the Time in Seconds for Execution of 
TRIGRD for a Case Having n Points 


n 

t 

t/Cn'^/S) 

t/(n log n) 

25 

0.061 

0.00084 

0.00175 

100 

0.346 

0.00075 

0.00173 

200 

0.951 

0.00081 

0.0020? 

500 

2.211 

0.00056 

0.00164 


The values at p of the first partial derivatives of this 
quadratic are stored as estimates of 9z/9x and 9z/9y at p. Execution 
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time on the Uni vac 1108 averages 8 milliseconds per point at which 
partlals are estimated. 

This method of estimating derivatives is the most ad hoc 
part of our entire surface interpolation method. We intend to investigate 
the effect of various parametric and algorithmic changes in this procedure. 


2.H LOOKUP IN THE TRIANGULAR GRID 

Given an arbitrary point q = (x, y) and an index k in the 
range 1 i k i n^i where n^ is the total number of triangles, the subroutine 
TRFINO tests to see if q is in the triangle whose index is k. 

If so, the index k is returned. If not, q must be outside 
one of the edges of the triangle. In this case, TRFIND resets k to be the 
index of the neighboring triangle on the other side of that edge and loops 
back to test q in this new triangle k. If there is no neighboring 
triangle, the fact that q is outside the triangular grid is reported. 

This approach is particularly efficient for the case of 
interpolating to points of a rectangular grid, since t!:e search can always 
be started at the triangle in which the previous point was found. When 
a new row of the rectangular grid is started, the search can be started 
in the triangle in which the first point of the previous row was found. 


2.5 INTERPOLATION IN A TRIANGLE 

The interpolation subroutine TVAL makes use of the piecewise 
cubic macroelement of Clough and Tocher (1965). A tutorial derivation 
of this element and a discussion of some alternative ways to organize 
its computation are given in Lawson (1976a). Quadrature properties 
of the element are derived in Lawson (1976b). 

Definition of this element involves partitioning the triangle 
into three subtriangles by drawing internal boundaries from the centroid 
to each vertex. In each of these three subtriangles the element is 
a cubic polynomial in x and y 


3 3-i 

z = E E aij xi yj 
i=0 J=0 


The element matches nine items of data, the function value 
and first partials with respect to x and y at the three vertices. 

It has continuity across the internal and external boundaries of 
the triangle. It is exact for quadratic data. 

Since it is a piecewise cubic, it is straightforward to ob- 
tain expressions for computing its first partial derivatives. The sub- 
routine TVAL includes an option to compute dz/dx and dz/dy as well as z 
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Starting with the inforination x, y, z, dz/dx, and 8z/8y 
at the verticea of a triangle, this method requires 55 multiplies, 

65 adds, and 4 divides to compute one Interpolated value. There are 
various possibilities for saving computed quantities that depend only 
on the triangle's data in order to out down the time to interpolate 
for a number of points in the same triangle. 

Execution time on the Uni vac 1108 averages 750 microseconds 
per interpolated point. This is about 10 times the cost of EXP or SIN. 

This timing of the interpolation subroutine TVAL Includes 
execution of the look-up subroutine TRFIND, usually beginning the look- 
up in the correct triangle or an immediately neighboring triangle. 


2.6 EXAMPLES 

Figure 2-3 shows a set S consisting of 26 points in the 
plane. Figure 2-4 is the triangular grid constructed for the set S 
by TRIGRD. In view of the results of Section 3.2, this is a Thlessen 
triangulation for S. 
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Figure 2-3* Set of 26 Points for Examples 1 and 2 
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Figure 2>4. Thiessen Triangular Grid for Examples 1 and 2 


2.6.1 Quadratic Test Case 

Values are assigned at the points of Figure 2-3 by computing 
2 - (-1 + 2x - 3y ♦ 4x^ ~ xy * 9y^)/8. Using ESTPD to estimate first 
partial derivatives and TRFIND and TVAL to interpolate to points of 
a 51 X 51 point rectangular grid, we then obtained the contour plot 
of Figure 2-5 and the perspective plot of Figure 2-6. This illustrates 
the exactness of the method for quadratic data. 

This case required 1.9 sec of Uni vac 1108 CPU time to build 
the triangular grid and Interpolate to the rectangular grid. It then 
used 15.3 sec in the plotting subroutines. 
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Figure 2-6. Perspective Plot for Example 1 
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Exponential Test Case 

For this case the z data is computed as 

2 2 , 

z s e“2(x ♦ y ) 


Again estimating partial derivatives at the 26 data points 
and then interpolating to a 51 x 51 rectangular grid, Figures 2-7 and 
2-8 are obtained. The TOst noticeable defect in the surface produced 
is the kink in the contour plot near (x,y) s (0.2, -0.4). This also 
appears as a groove in the perspective plot. From Figure 2-3, however, 
it is seen that this corresponds to a region in which there is a lack 
of data. Computer time was similar to the first test case. 



Figure 2-7. Contour Plot for Example 2 
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Figure 2-8. Perspective Plot for Example 2 


2.6.3 A Case With More Points 

The third test case is the same exponential function on 
a set of 500 points. The grid produced for this case has 985 triangles. 

This data was interpolated to a 21 x 21 point rectangi‘?'r 
grid for plotting (see Figures 2-9« 2-10, and 2-11). This case 
6.25 sec of CPU time to triangulate and interpolate. It used M.Oi 
sec in the plotting subroutines. The plotting was faster because of 
the much coarser rectangular grid. 
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Figure 2-11. Perspective Plot for Example 3 
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SECTION 111 

THEOHETICAL HESULTS ON Tftl ANGULATION 


3.1 THREb CRITERIA FOR TRIANGULATION OF A STRICTLI CONVEX 

gUALRlLATbHAL 

Ne will call a quadrllataral Q strictly convex if each 
of its four interior angles measures less than IbO^. Such a quadrilateral 
can be partitioned into two triangles in two possible ways. Three 
criteria will be described for choosing a preferred triangulation of 
a strictly convex quadrilateral. 


3.1.1 The Hax*Min Angle Criterion 

Choose the triangulation of Q that maximizes the minimum 
interior angle of the two resulting triangles. Either choice can be 
made in the case of a tie. For example, ir Figure 3-1 ^cab is the 
smallest angle in triangles fi and g^t x.cdb is the smallest angle in 
triangles f2 and g2, and ^edb is larger than ^tcab. Thus, the triangulation 
(b) in preferred over (a). 


3.1.2 The Circle Criterion 

Let K denote a circle passing through three of the vertices 
of a strictly convex quadrilateral Q. If the fourth vortex is Interior 
to K, insert the diagonal from this fourth vertex to the opposite vertex. 
If the fourth vertex is exterior to K, Insert the other diagonal. 

If the fourth vertex is on K, insert either diagonal (see Figure 3-2 
for an example). Note that when all four vertices are not on a common 
circle, the same triangulation will be selected regardless of which 
set of three vertices is used to construct the circle. 


3.1.3 The Thiessen Region Criterion 

Let Rg denote the closure of the region of the plane censistina 
of all points that are closer to point a than to points b^ ;, or d. 
Similarly define regions h^, R^, and R^j surrounding points b, c, and d, 
respectively. These regions are called 'Ihieasen regions following 
Powell (1976) and Rhynsburger (1973). Ihese proximity regions are 
also identified by other names in the mathematical literature. 

Two of the points a, b, c, or d will be called Thiessen 
neighbors if their Thiessen regions are in contact. They are strong 
Thiessen neighbors if the contact is along a line segment cf nonzero 
length. They are weak Thiessen neighbors if the contact i.s at cne 
point only. 
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Figure 3-1. The Max-Min Angle Criterion 



To triangulate a strictly convex quadrilateral Q, insert 
the oiagonal that connects a pair of strong Thiessen neighbors. A 
strictly convex quadrilateral can have at nost one pair of opposite 
vertices that are strong Thiessen neighbors. If neither pair of opposite 
vertices are strong Thiessen neighbors, then both pairs will be weak 
neighbors and either diagonal can be inserted (see Figure 3-3 for an 
example ) . 
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Figure 3-3* The Thiessen Region Criterion 


3.1.4 Equivalence of Thetc Three Criteria for Strictly Convex 

Quadrilaterals 

The first observation to be made about these three criteria 
is that they give identical results for strictly convex quadrilaterals. 
This can be verified by noting that all three criteria have the same 
neutral case and then studying perturbations from the neutral case. 

The neutral case for all three criteria is the case in 
which all four vertices of the quadrilateral lie on a common circle. 

To verify this last statement consider a quadrilateral 
Q whose vertices a, b, c, and d all lie on a common circle K (see Figure 
3-4). Suppose arc be is shorter than arcs cd, da, or ab. If the angular 
measure of arc be is 26, then angles cab and cdb are each of measure 
0, and these two angles are each the minimum angle for one of the possible 
triangulations. Thus, this is a tie case for the max-min angle criterion. 



Figure 3-4. The Neutral Case for the Max-Min Angle Criterion 
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Constructing Thiessen regions for the case of four points 
on a common circle results in the four Thiessen regions meeting at the 
center of the circle, as in Figure 3-5. Thus, each pair of opposite 
vertices are weak Thiessen neighbors, which is the neutral care for 
this criterion. 

Further analysis, the details of which will be omitted, 
shows that moving one point, say point d in Figures 3-^ and 3-5, inside 
circle K causes ^ cdb to become larger and causes points b and d to 
become strong Thiessen _neighbors. Thus, all three criteria will choose 
to introduce the edge bd. 


3.1.5 A Local Optimization Procedure 

Let e denote an Internal edge in a triangular grid T. 

Consider the quadrilateral Q formed by the two triangles having e as a 
common edge. If Q is not strictly convex then e cannot be considered for 
swapping. Otherwise, if Q is strictly convex, apply any one of the three 
equivalent criteria discussed in the preceding sections. Replace e by 
the other diagonal of Q if this is preferred by the criterion. Otherwise 
if e is preferred or if the decision is neutral, leave e as it is. 

If the criterion used is the circle test and if the circle 
used is the clrcumcircle of one of the triangles containing the edge 
e, then it is not necessary to do an initial test for Q being strictly 
convex since the correct decision will be made anyway. That is, if 
Q is not strictly convex, then the clrcumcircle of one triangle will 
not enclose the vertex of the other triangle so the decision will be 
not to swap the edge e. 

An internal edge of a triangulation T will be called locally 
optimal if application of the local optimization procedure to it would 
not swap it. 



Figure 3-5. The Neutral Cjise for the Thiessen Region Criterion 
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3.2 GLOBAL CONSEQUENCES OF THE LOCAL OPTIMIZATION PROCEDURE 

The local optimization procedure has a number of consequences 
that are useful in suggesting, and proving properties of, a variety 
of possible triangulation algorithms. 


3.2.1 A Linear Ordering of Triangulations 

Let S be a set of n points in the plane and let ^ denote 
the set of all triangulations of S. As has previously been noted, 
all triangulations Tc^ have the same number of triangles, say n^. 

With each Tc.^ associate an indicator vector of n^ components 
constructed as follows: Determine the measure of the smallest angle 

in each of the n^ triangles of T and sort these angular measures in 
nondecreasing order. The triangulations in ^9" can then be linearly 
ordered by the lexicographic ordering of their associated vectors. 


3. 2. 1.1 Iheorea. Let T be a trlangulatlon of a finite point set 

S. LfiL e be an Internal edge of T. Suppose aPDlleatlon of the local 
optiliaation procedure to e leads to a -swap, replacing e bv a new edge 
e', and thus producing a new triangulation t* qL s. Then t < t’; i.e. . 

T' strictly J'ollows T in the linear ordering defined above. 

Proof . Let v be the indicator vector for T. The measures 
of the smallest angles in the two triangles in T sharing the edge e 
occur as two of the components of v, say Vj and vjj, with j < k and 
thus Vj 1 vt(. Since a swap was made when e was tested, the smallest 
angles in both of the new triangles of T' sharing the edge e' must 
be strictly greater than vj. It follows that the indicator vector 
v» for T' strictly follows v in lexicographic order and thus T < T’.l 

This theorem can be used to prove finite termination for 
a variety of possible triangulation algorithms that repeatedly apply 
the local optimization procedure to all internal edges of a sequence 
of triangulat.''ons. Since there are only a finite number of possible 
triangulations of a point set S and each swap produced by the local 
optimization procedure causes a strict advance through a linear ordering 
of the triangulations, it follows that after some finite number of 
swaps a triangulation T* will be reached such that each internal edge 
in T* will be left unswapped when tested by the local optimization 
procedure; i.e., all internal edges are locally optimal. 


3. 2;. 1.2 Theorem. All internal edges of a trlangulatlon T of a 

finite point set s are locally optimal if and only if do point of s 
is interior to any circumcircle _of a triangle of T. 

Proof of "if' . Assume no point of S is interior to any 
circumcircle of a triangle of T. Consider the application of the local 
optimization procedure to any internal edge e in T. Let f and g denote 
the two triangles sharing the edge e. Let d denote the vertex of g 
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opposite to edge e. By hypothesis d Is not interior to the circumcircle 
of triangle f. Thus the local optimization procedure will not swap 
e; that is, e is already locally optimal. 

Proof of "only if" . Assume all internal edges of T are 
locally optimal. Suppose the theorem is false; i.e., suppose there 
is a triangle f in T with vertices a, b, and c and circumcircle K such 
that there is a point p of S interior to K. 

Without loss of generality assume edge ac is the nearest 
edge of^abc to p as in Figure 3-6. Denote the perpendicular distance 
from ac to p by 6 . Among all triangles of T, whose circumcircle contains 
p as an interior point, assume without loss of generality that none 
is at a distance of less than 6 from p. 

Since p is on the opposite side of ac ^rom b, the edge 
ac is not a boundary edge of T. Thus there is another triangle, say 
^cq, sharing the edge ac. The vertex q cannot be interior to the 
circle K as this would contradict the hypothesis that edge ac is locally 
optimal. Thus q is on K or exterior to K. 

Without loss of generality, assume edge cq is the closest 
edge of Aacq to p as in Figure 3-7. Note that the distance from cq 
to p is less than 6. Thus a contradiction will be reached if it is 
shown that the circumcircle K* of Aacq contains p as an interior point. 
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q 



Figure 3-7. Theorem 3. 2. 1.1 


If q is on K, then K' = K, and so p lies interior to K'. 

If q is outside K, then K' intersects K only at a and c and encloses 
all of the interior of K to the right of ac in Figure 3-7. In particular 
K* encloses p. ■ 


3.2.2 The Global Thiessen Triangulation 

The notion of Thiessen regions and an associated triangulation 
v.’erfc defined for strictly convex quadrilaterals in Section 3.1.3. 

The same definitions can be applied to any finite set of points in 
the plane. First the Thiessen region surrounding each point can be 
determined. Then line segments can be drawn connecting points that 
are strong Thiessen neighbors. This will provide a polygonal grid. 

For points in general position, the grid will in fact consist 
of triangles. Any polygon of the grid that is not a triangle will 
have all of its vertices on a circle and all nonadjacent pairs of vertices 
will be weak Thiessen neighbors. Any such k-point polygon can be triangulated 
by connecting any k-3 pairs of its vertices as long as no crossing 
lines are drawn. A triangulation in which all strong Thiessen neighbors 
are connected will be called a Thiessen trlangulatlon . 


£.2.2.1 Lemma . Let S be a set p^f n points in the plane and let 
S' be a subset of s. Two points of S that are Thiessen neighbors in 
S remain so in s and if the y are st rong Thiess en neighbors in S Uiey. 
remain strong Thiessen neighbors in ?. 

Proof . Consider the effect of removing one point, say 
p, from S, leaving a subset S'. The only change that takes place in 
transforming the Thiessen regions for S to form the Thiessen regions 
for S' is a redistribution of the part of the plane that was the Thiessen 
region for p. This region will be partitioned, with various portions 
being absorbed into the neighboring Thiessen regions. In this process 
no boundaries between pairs of remaining Thiessen regions are shortened. 
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Thus neighbors remain neighbors and strong neighbors remain strong 
neighbors. Clearly the same is true for the removal of any number 
of points from a finite set, as the removals can be done one at a time. I 


3. 2 . 2. 2 Ibaorea. All internal edges of a triangulation T o£_a 

finite point set S are locally optimal if and only if T ia a Ihieasen 

trlaMulatlon of s. 

Proof of *‘if»' . Assume T is a Thiessen triangulation of S. 

Let e be an Internal edge of T connecting vertices a and c and belonging 
to triangles Aabc andAcda. If the quadrilateral Q s abed is not strictly 
convex, then e cannot be swapped and is thus locally optimal. 

Consider then the case of Q being strictly convex. By 
hypothesis, a and c are Thiessen neighbors in S. By Lemma 3.2.2. 1 they 
are also Thiessen neighbors relative to the point set |a, b, c, d|. With 
the quadrilateral Q = abed being strictly convex, this implies e is 
locally optimal. 

Proof of "only if*' . Assume all internal edges of T are 
locally optimal. Suppose the theorem is false. Then there is some 
pair of strong Thiessen neighbors in S, say points p and q, that are 
not connected by an edge in the triangulation T. 

Define B to be the polygonal curve whose constituent line 
segments are the segments that occur as edges opposite vertex q in 
those triangles that have q as a vertex. If q is not a boundary point 
of T, then B is a closed polygon with q in its interior and p lying 
exterior to it. Clearly a line segment from p to q would Intersect B. 

If q is a boundary point of T, then B is an open polygonal 
curve with end points on the boundary of T at the two points immediately 
adjacent to q on each side of q along the boundary. Although in this 
case B does not surround q, it still follows that ^ must intersect 
P owing to the convexity of the region covered by T. 

Since p and q are strong Thiessen neighbors in S, there 
can be no other points of S on the line segment Thus the Intersection 
of ^ with B is not at a point of S on B but must be strictly between 
a pair of points of S on B, say points r and s. 

Thus the triangle Aqrs is a triangle of T having the property 
that r and s are on strictly opposite sides of ^ and p and q are on 
strictly opposite sides of rs, as is illustrated in Figure 3-8. By hypo- 
thesis, all internal edges of T are locally optimal. By Theorem 3. 2. 1.2 
this implies that point p is not in the interior of the clrcumcircle K of 
Aqrs. From the equivalence of the circle test and the Thiessen criterion 
for strictly convex quadrilaterals, this would imply that p is not a 
strong Thiessen neighbor of q relative to the point set | p, s, q, r|. 

By Lemma 3.2.2. 1, however, p and q are strong Thiessen 
neighbors relative to |p, s, q, r| since they are strong Thiessen neigh- 
bors in S. Thus a contradiction has been reached. I 
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Figure 3-8. Theorem 3. 1.2. 2 


3.3 McLAlN'S TRIANGULATION METHOD 

The triangulation method described in McLain (1976) builds 
its grid one triangle at a time in such a way that each triangle con- 
structed is a triangle of the final grid. This is in contrast to methods 
that involve triangle modification steps such as are used in Section 
2.2, in Lawson (1972), and by frank Little of the University of Utah 
CAGD group. 


The paper, McLain (1976), with the subsequent errata, leaves 
open the question of the characterization of the grid the algorithm 
produces. We find that the results of the preceding sections can be 
used to show that the grid produced by McLain's method is in fact a 
Thiessen grid. 

Let S be a set of n points. Define Tq to be a single edge 
belonging to some Thiessen triangulation for S. For example, 1 q could 
be a boundary edge of the convex hull of S. For k 2 1 define 
be a configuration of k triangles that is a subset of some Thiessen 
triangulation of S and contains T^-i as a subset. In general, the 
configurations T|^ are not necessarily convex. 

Given some Tij, how can one more triangle of a Thiessen 
triangulation of S be found to advance to 

Let ab be an edge belonging to just one triangle, say ^abc, 
in Tjj. Let be the _subset of S consisting of the points lying on 
the oppo^te side of ab from c. (If there are none, then try another 
edge as ab or terminate.) 

From our inductive assumption that Ti< is a subset of a 
Thiessen trianguiation of .S, there must be a point p in such that 
adjoining Aabp to Tj^ gives a configuration T^^+i that is also a subset 
of a Thiessen trianguiation of S. 

By Theorems 3.2. 1.2 and 3. 2. 2. 2 we know that such a triangle, 
Aabp, must not contain any points of S|^ in its interior, ouch a triangle 
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is Just what is selected by Moldin' s method, since he selects p such 
that the signed distance from ab of the center of the circumclrcle 
of Aabp is the algebraically smallest possible among all choices of 
p in S|(. The signed distance from aB is positive on the side of aB 
opposite to c. 


3.^ LIMITS ON GRID CHANGES WHEN ADDING A NEW POINT 

This section presents results that limit the amount of 
edge testing needed in algorithms such as ours in Section 2.2 that 
transform from the Thiessen triangulation of one point set to that 
for the set augmented by one new point. 

Let S|).^ be a set of n-1 points in the plane. Let p be 
a point not in Sn-i and define S„ - Sn-i U |pf* ^at Tn-1 Pa a Thiessen 
triangulation for S„_^ . 

Given Tn- 1 , initial triangulation Tn^^^ for Sn can be 
constructed by inserting all edges that connect the new point p to 
points of Sn -1 without crossing edges already present in Tn-i- A special 
case arises if p falls on an edge, say ac, already present in T^_i . 

Then the edge ac must be replaced by the two edges ap and For 

our theoretical discussion it will be easier to assume p does not fall 
on an edge of Tn-i. The case of p arbitrarily close to an edge of 
Tn _1 is of course permitted, and analysis of this case can be used 
to justify the replacement of ac by ap and 

An edge e in a triangulation will be called converged if 
it is locally optimal in the present trlangulatlon and if in addition 
it can be proved that no sequence of applications of the local optimiza- 
tion procedure to the various edges could lead to a decision to swap e. 

Assuming p does not lie on any edge of Tj^.i, it will be 
shown that all of the initial edges inserted connecting p to points 
of Sn -1 are converged. Any edge opposite to p in a triangle must be 
tested once using the local optimization procedure. If it remains 
unchanged after testing, then it is converged. If it is swapped by 
the procedure, then the new edge introduced is converged and the two 
edges opposite p in the two new triangles must each be tested. 


3.^.1 Theorem 

Aaauae P is strictly outside the convex hull of Sn_i aoil 
Tn^^^ la formed by goppecting P to all boundary points of T ^-1 LlaL 
can be reached without cross ing any edges of T^_i . Then all of the 

aew- edgea . are converged . 

Proof . Let ^ be a new edge connecting p to a point q 
on the boundary of Tn-i. If in the course of triangulating Sr It is 
ever to be decided to swap ^ for some other edge ab, then a and b 
must (at least) be points of Sn-i such that a and b are on strictly 
opposite sides of pq and p and q are on strictly opposite sides of ab. 
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This is impossible since the line segment ab could not 
pass outside the convex hull of and thus could not intersect ^ 
strictly between the boundary point q and the ooint p which is exterior 
to the convex hull of Sp.i. I 


3.^.2 Theorem 

Aaaume P is interior to the convex hull of Sn.i and .In 
fact iaterlor tfl..agae -triangle Aabc n£ Tq-i. Assume Tn(l) is formed 
t?y connecting P a, b, Mil O. These three new edges are c onverged . 

Pranf « since T*_i is a Thiessen triangulation of S^-i, 
the circumcircle K of Aabc contains no points of Sn.i in its interior. 

If in the course of triangulating Sjj^it is ever decided to swap 
for instance, for josie other edge, i^, then r and s must not be interior 
to K, r and s mur t be on strictly opposite sides of pa, p and a must 
be on strictly opposite sides of rs, and a must be strictly outside 
the circle K' through r, p, and s (see Figure 3.9). 

Partition circle K‘ into the three arcs rp, and sr. 

Note that arc fp intersects circle K since p is inside K and r is on 
or outside K. Call this intersection point r'. Similarly, arc ps 
intersects circle K, say at a point s'. Then arc r^^' of K* is inside 
K oecause p is inside K and a»*e intersection points of K* 

with K. It follows that arc s^sr^ of K' is outside K. The arc of 
K between s' and r' that lies Inside K' contains the point a. Thus 
it is impossible for a to be outside K' as it must be for pa to be 
swapped. I 


3 . 4.3 Theorem 

Let Tp(l) be a trianeulation of Sp = Sp-i U|p|» Lfit ^bp 
M^i^abc fcs adJasent triangles of T^d) and assume Aabc was also a 
triangle of Tj^.^. Suppose apDllcatlon^of the_looal optimiz^ion procedure 
to edge He leads to a swan. r.eplaclng be Jiy oa. Then edge pa la CQnvgCgfi4« 

Proof . Note that the symbols, a, b, and e do not necessarily 
denote the first vertices to which p was connected as was tlie case 
in Theorem 3.4.2. The notation for this theorem was selected, however, to 
permit the proof to be identical to that of Theorem 3*4.2 (see Figure 3.10).i 


3-11 


1 


I 


x( „ U I I 1 I 1 

77-30 


c 



Figure 3-9. Theorem 3. **.2 



Figure 3-10. Theorem 3- *1.3 


3.5 CONCLUSIONS 

Triangular grid construction is certainly not as well understood 
as sorting of scalars, but at least a general notion of what to expect 
in running time, storage usage, and properties of the final triangulation 
is now available. Thus an algorithm with an 0(n2) time estimate must 
be regarded as inefficient except possibly for small n. There remain 
a wide variety of possible triangulation algorithms with time estimates 
in the neighborhood of 0(n^'^3), it will require more time, experience 
and direct comparative testing to sort these out. 

C^ interpolation on triangles is still very ad hoc. Three 
different methods are used by Akima ( 1975), McLain (19^), and the 
present paper. It appears that our method requires the least amount 
of auxiliary stored information per node (two first partial derivatives 
per node) and the fewest number of operations per interpolated value 
once the auxiliary nodal information has been computed and stored. 


3-12 


77-30 


Direct comparative testa would be needed, however, to assess accuracy 
and actual execution times. 

There is much scope for additional work on this problem, 
including generalizations such as smoothing instead of interpolation, 
the introductions of constraints, and permitting the domain of the 
independent variables to be the surface of a sphere or three-dimensional 
space instead of the plane. 
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